<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>tanh_sinh</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Math Toolkit 3.3.0">
<link rel="up" href="../double_exponential.html" title="Double-exponential quadrature">
<link rel="prev" href="de_overview.html" title="Overview">
<link rel="next" href="de_tanh_sinh_2_arg.html" title="Handling functions with large features near an endpoint with tanh-sinh quadrature">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="de_overview.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="de_tanh_sinh_2_arg.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.double_exponential.de_tanh_sinh"></a><a class="link" href="de_tanh_sinh.html" title="tanh_sinh">tanh_sinh</a>
</h3></div></div></div>
<pre class="programlisting"><span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">&gt;</span>
<span class="keyword">class</span> <span class="identifier">tanh_sinh</span>
<span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
    <span class="identifier">tanh_sinh</span><span class="special">(</span><span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">Real</span><span class="special">&amp;</span> <span class="identifier">min_complement</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">min_value</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()</span> <span class="special">*</span> <span class="number">4</span><span class="special">)</span>

    <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
    <span class="keyword">auto</span> <span class="identifier">integrate</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span>
                   <span class="identifier">Real</span> <span class="identifier">tolerance</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
                   <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
                   <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
                   <span class="identifier">std</span><span class="special">::</span><span class="identifier">size_t</span><span class="special">*</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()))</span> <span class="keyword">const</span><span class="special">;</span>

    <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
    <span class="keyword">auto</span> <span class="identifier">integrate</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span>
                   <span class="identifier">tolerance</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
                   <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
                   <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
                   <span class="identifier">std</span><span class="special">::</span><span class="identifier">size_t</span><span class="special">*</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()))</span> <span class="keyword">const</span><span class="special">;</span>

<span class="special">};</span>
</pre>
<p>
        The <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature routine provided by boost
        is a rapidly convergent numerical integration scheme for holomorphic integrands.
        By this we mean that the integrand is the restriction to the real line of
        a complex-differentiable function which is bounded on the interior of the
        unit disk <span class="emphasis"><em>|z| &lt; 1</em></span>, so that it lies within the so-called
        <a href="https://en.wikipedia.org/wiki/Hardy_space" target="_top">Hardy space</a>.
        If your integrand obeys these conditions, it can be shown that <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code>
        integration is optimal, in the sense that it requires the fewest function
        evaluations for a given accuracy of any quadrature algorithm for a random
        element from the Hardy space.
      </p>
<p>
        A basic example of how to use the <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature
        is shown below:
      </p>
<pre class="programlisting"><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">5</span><span class="special">*</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">7</span><span class="special">;</span> <span class="special">};</span>
<span class="comment">// Integrate over native bounds of (-1,1):</span>
<span class="keyword">double</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">);</span>
<span class="comment">// Integrate over (0,1.1) instead:</span>
<span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="number">1.1</span><span class="special">);</span>
</pre>
<p>
        The basic idea of <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature is that a variable transformation
        can cause the endpoint derivatives to decay rapidly. When the derivatives
        at the endpoints decay much faster than the Bernoulli numbers grow, the Euler-Maclaurin
        summation formula tells us that simple trapezoidal quadrature converges faster
        than any power of <span class="emphasis"><em>h</em></span>. That means the number of correct
        digits of the result should roughly double with each new level of integration
        (halving of <span class="emphasis"><em>h</em></span>), Hence the default termination condition
        for integration is usually set to the square root of machine epsilon. Most
        well-behaved integrals should converge to full machine precision with this
        termination condition, and in 6 or fewer levels at double precision, or 7
        or fewer levels for quad precision.
      </p>
<p>
        One very nice property of tanh-sinh quadrature is that it can handle singularities
        at the endpoints of the integration domain. For instance, the following integrand,
        singular at both endpoints, can be efficiently evaluated to 100 binary digits:
      </p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)*</span><span class="identifier">log1p</span><span class="special">(-</span><span class="identifier">x</span><span class="special">);</span> <span class="special">};</span>
<span class="identifier">Real</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="special">(</span><span class="identifier">Real</span><span class="special">)</span> <span class="number">0</span><span class="special">,</span> <span class="special">(</span><span class="identifier">Real</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
</pre>
<p>
        Now onto the caveats: As stated before, the integrands must lie in a Hardy
        space to ensure rapid convergence. Attempting to integrate a function which
        is not bounded on the unit disk by tanh-sinh can lead to very slow convergence.
        For example, take the Runge function:
      </p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f1</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">1</span><span class="special">+</span><span class="number">25</span><span class="special">*</span><span class="identifier">t</span><span class="special">*</span><span class="identifier">t</span><span class="special">);</span> <span class="special">};</span>
<span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="special">-</span><span class="number">1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
</pre>
<p>
        This function has poles at ± ⅈ/5, and as such it is not bounded
        on the unit disk. However, the related function
      </p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f2</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">1</span><span class="special">+</span><span class="number">0.04</span><span class="special">*</span><span class="identifier">t</span><span class="special">*</span><span class="identifier">t</span><span class="special">);</span> <span class="special">};</span>
<span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f2</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="special">-</span><span class="number">1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
</pre>
<p>
        has poles outside the unit disk (at ± 5ⅈ), and is therefore in
        the Hardy space. Our benchmarks show that the second integration is performed
        22x faster than the first! If you do not understand the structure of your
        integrand in the complex plane, do performance testing before deployment.
      </p>
<p>
        Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate
        of the L<sub>1</sub> norm of the integral along with the requested integral. This is
        to establish a scale against which to measure the tolerance, and to provide
        an estimate of the condition number of the summation. This can be queried
        as follows:
      </p>
<pre class="programlisting"><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">5</span><span class="special">*</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">7</span><span class="special">;</span> <span class="special">};</span>
<span class="keyword">double</span> <span class="identifier">termination</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">());</span>
<span class="keyword">double</span> <span class="identifier">error</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</span>
<span class="identifier">size_t</span> <span class="identifier">levels</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="number">1.0</span><span class="special">,</span> <span class="identifier">termination</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">L1</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">levels</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">condition_number</span> <span class="special">=</span> <span class="identifier">L1</span><span class="special">/</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">abs</span><span class="special">(</span><span class="identifier">Q</span><span class="special">);</span>
</pre>
<p>
        If the condition number is large, the computed integral is worthless: typically
        one can assume that Q has lost one digit of precision when the condition
        number of O(10^Q). The returned error term is not the actual error in the
        result, but merely an a posteriori error estimate. It is the absolute difference
        between the last two approximations, and for well behaved integrals, the
        actual error should be very much smaller than this. The following table illustrates
        how the errors and conditioning vary for few sample integrals, in each case
        the termination condition was set to the square root of epsilon, and all
        tests were conducted in double precision:
      </p>
<div class="informaltable"><table class="table">
<colgroup>
<col>
<col>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
                <p>
                  Integral
                </p>
              </th>
<th>
                <p>
                  Range
                </p>
              </th>
<th>
                <p>
                  Error
                </p>
              </th>
<th>
                <p>
                  Actual measured error
                </p>
              </th>
<th>
                <p>
                  Levels
                </p>
              </th>
<th>
                <p>
                  Condition Number
                </p>
              </th>
<th>
                <p>
                  Comments
                </p>
              </th>
</tr></thead>
<tbody>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="number">5</span> <span class="special">*</span>
                  <span class="identifier">x</span> <span class="special">+</span>
                  <span class="number">7</span></code>
                </p>
              </td>
<td>
                <p>
                  (0,1)
                </p>
              </td>
<td>
                <p>
                  3.5e-15
                </p>
              </td>
<td>
                <p>
                  0
                </p>
              </td>
<td>
                <p>
                  5
                </p>
              </td>
<td>
                <p>
                  1
                </p>
              </td>
<td>
                <p>
                  This trivial case shows just how accurate these methods can be.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
                  <span class="special">*</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
                </p>
              </td>
<td>
                <p>
                  0, 1)
                </p>
              </td>
<td>
                <p>
                  0
                </p>
              </td>
<td>
                <p>
                  0
                </p>
              </td>
<td>
                <p>
                  5
                </p>
              </td>
<td>
                <p>
                  1
                </p>
              </td>
<td>
                <p>
                  This is an example of an integral that Gaussian integrators fail
                  to handle.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span>
                  <span class="special">/</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
                </p>
              </td>
<td>
                <p>
                  (0,+∞)
                </p>
              </td>
<td>
                <p>
                  8.0e-10
                </p>
              </td>
<td>
                <p>
                  1.1e-15
                </p>
              </td>
<td>
                <p>
                  5
                </p>
              </td>
<td>
                <p>
                  1
                </p>
              </td>
<td>
                <p>
                  Gaussian integrators typically fail to handle the singularities
                  at the endpoints of this one.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">x</span> <span class="special">*</span>
                  <span class="identifier">sin</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">exp</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">exp</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">))))</span></code>
                </p>
              </td>
<td>
                <p>
                  (-1,1)
                </p>
              </td>
<td>
                <p>
                  7.2e-16
                </p>
              </td>
<td>
                <p>
                  4.9e-17
                </p>
              </td>
<td>
                <p>
                  9
                </p>
              </td>
<td>
                <p>
                  1.89
                </p>
              </td>
<td>
                <p>
                  This is a truly horrible integral that oscillates wildly and unpredictably
                  with some very sharp "spikes" in it's graph. The higher
                  number of levels used reflects the difficulty of sampling the more
                  extreme features.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span>
                  <span class="number">0</span> <span class="special">?</span>
                  <span class="number">1</span> <span class="special">:</span>
                  <span class="identifier">sin</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
                  <span class="special">/</span> <span class="identifier">x</span></code>
                </p>
              </td>
<td>
                <p>
                  (-∞, ∞)
                </p>
              </td>
<td>
                <p>
                  3.0e-1
                </p>
              </td>
<td>
                <p>
                  4.0e-1
                </p>
              </td>
<td>
                <p>
                  15
                </p>
              </td>
<td>
                <p>
                  159
                </p>
              </td>
<td>
                <p>
                  This highly oscillatory integral isn't handled at all well by tanh-sinh
                  quadrature: there is so much cancellation in the sum that the result
                  is essentially worthless. The argument transformation of the infinite
                  integral behaves somewhat badly as well, in fact we do <span class="emphasis"><em>slightly</em></span>
                  better integrating over 2 symmetrical and large finite limits.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span>
                  <span class="special">(</span><span class="number">1</span>
                  <span class="special">-</span> <span class="identifier">x</span>
                  <span class="special">*</span> <span class="identifier">x</span><span class="special">))</span></code>
                </p>
              </td>
<td>
                <p>
                  (0,1)
                </p>
              </td>
<td>
                <p>
                  1e-8
                </p>
              </td>
<td>
                <p>
                  1e-8
                </p>
              </td>
<td>
                <p>
                  5
                </p>
              </td>
<td>
                <p>
                  1
                </p>
              </td>
<td>
                <p>
                  This an example of an integral that has all its area close to a
                  non-zero endpoint, the problem here is that the function being
                  integrated returns "garbage" values for x very close
                  to 1. We can easily fix this issue by passing a 2 argument functor
                  to the integrator: the second argument gives the distance to the
                  nearest endpoint, and we can use that information to return accurate
                  values, and thus fix the integral calculation.
                </p>
              </td>
</tr>
<tr>
<td>
                <p>
                  <code class="computeroutput"><span class="identifier">x</span> <span class="special">&lt;</span>
                  <span class="number">0.5</span> <span class="special">?</span>
                  <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
                  <span class="special">/</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">x</span>
                  <span class="special">*</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span>
                  <span class="special">((</span><span class="identifier">x</span>
                  <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">xc</span><span class="special">)))</span></code>
                </p>
              </td>
<td>
                <p>
                  (0,1)
                </p>
              </td>
<td>
                <p>
                  0
                </p>
              </td>
<td>
                <p>
                  0
                </p>
              </td>
<td>
                <p>
                  5
                </p>
              </td>
<td>
                <p>
                  1
                </p>
              </td>
<td>
                <p>
                  This is the 2-argument version of the previous integral, the second
                  argument <span class="emphasis"><em>xc</em></span> is <code class="computeroutput"><span class="number">1</span><span class="special">-</span><span class="identifier">x</span></code>
                  in this case, and we use 1-x<sup>2</sup> == (1-x)(1+x) to calculate 1-x<sup>2</sup> with
                  greater accuracy.
                </p>
              </td>
</tr>
</tbody>
</table></div>
<p>
        Although the <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature can compute integral over
        infinite domains by variable transformations, these transformations can create
        a very poorly behaved integrand. For this reason, double-exponential variable
        transformations have been provided that allow stable computation over infinite
        domains; these being the exp-sinh and sinh-sinh quadrature.
      </p>
<h5>
<a name="math_toolkit.double_exponential.de_tanh_sinh.h0"></a>
        <span class="phrase"><a name="math_toolkit.double_exponential.de_tanh_sinh.complex_integrals"></a></span><a class="link" href="de_tanh_sinh.html#math_toolkit.double_exponential.de_tanh_sinh.complex_integrals">Complex
        integrals</a>
      </h5>
<p>
        The <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> integrator
        supports integration of functions which return complex results, for example
        the sine-integral <code class="computeroutput"><span class="identifier">Si</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span></code> has
        the integral representation:
      </p>
<div class="blockquote"><blockquote class="blockquote"><p>
          <span class="inlinemediaobject"><img src="../../../equations/sine_integral.svg"></span>

        </p></blockquote></div>
<p>
        Which we can code up directly as:
      </p>
<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">Complex</span><span class="special">&gt;</span>
<span class="identifier">Complex</span> <span class="identifier">Si</span><span class="special">(</span><span class="identifier">Complex</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
   <span class="keyword">typedef</span> <span class="keyword">typename</span> <span class="identifier">Complex</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">value_type</span><span class="special">;</span>
   <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">;</span>  <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special">;</span>
   <span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;</span><span class="identifier">z</span><span class="special">](</span><span class="identifier">value_type</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="special">-</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">z</span> <span class="special">*</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">t</span><span class="special">))</span> <span class="special">*</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">z</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span><span class="identifier">t</span><span class="special">));</span> <span class="special">};</span>
   <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">quadrature</span><span class="special">::</span><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="identifier">value_type</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
   <span class="keyword">return</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">half_pi</span><span class="special">&lt;</span><span class="identifier">value_type</span><span class="special">&gt;())</span> <span class="special">+</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">half_pi</span><span class="special">&lt;</span><span class="identifier">value_type</span><span class="special">&gt;();</span>
<span class="special">}</span>
</pre>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
      Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
      Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
      Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
      Walker and Xiaogang Zhang<p>
        Distributed under the Boost Software License, Version 1.0. (See accompanying
        file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
      </p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="de_overview.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="de_tanh_sinh_2_arg.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
